function IR = LPM1(I)
 [F C] = size(I);
 
 Om = F/2;                  % coordenadas del centro de la imagen
 On = C/2;
 sx = F/2;                  % factores de escala
 sy = C/2;
 
 M=64;
 N=360;
  
 %img=zeros(M*2,N);
 
 R=0.4;
 a=2^(1/M);
 rmin=(a^(-M/2))*R;
 rmax=(a^(M/2))*R;

  cont=zeros(M*2,N);
  suma=zeros(M*2,N);
  img=zeros(M*2,N);
 
 
for xi = 1:F
  for yi = 1:C
    x = abs((xi - Om)/sx);
    y = abs(((-1)*yi + On)/sy);
     
   % x = (xi - Om)/sx;
   % y = (yi - On)/sy;
     
     r= sqrt(x^2 + y^2);

     if (r >=rmin) && (r <=rmax)
        
        theta = atan2(y,x);
        
        l1 =abs(floor((log2(r)-log2(R))/log2(a))+M/2);
        l2 =abs(floor((N*theta)/pi));
          
        if(l1>=0 && l1<2*M) && (l2>=0 && l2<N)
          img(l1+1,l2+1) = I(xi+1,yi+1);
           
          cont(l1+1,l2+1)=cont(l1+1,l2+1)+1;
          suma(l1+1,l2+1)=suma(l1+1,l2+1)+I(xi+1,yi+1);

        end
      end 
  end
end
cont(cont(:)==0)=1;
img(:)=suma(:)./cont(:);
IR=img;
end